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ABSTRACT 


An iterative numerical method- is studied for 
computing steady, three-dimensional, viscous, compres- 
sible flow fields about aerodynamic bodies . In order 
to develop the iterative method economically, the primary 
emphasis of this research was directed towards super- 
sonic, axisymmetric flow. However, the technique readily 
extends to three spatial dimensions. The viscous flow 
field about a cone-cylinder-flare body was calculated 
and compared to existing experimental data. Numerical 
predictions of the cone boundary layer and the flow 
field shock wave structure agreed with corresponding 
measurements. A separation was calculated at the 
cylinder-flare junction in six iterations; however, the 
size of the vortex did not correspond to the measured 
size. It was estimated that fifty iterations would be 
required to properly define' the vortex. Conclusions of 
the study are that the iteration technique is of limited 
value for plane two-dimensional and axisymmetric flows , 
but of great value for three-dimensional flows. 
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NOMENCLATURE 


Radius of the cylindrical cross-section 
Local sound speed 
Chapman-Rubesin constant 

Diameter of cylindrical portion of cone-cylinder- 
flare body 

Specific internal energy 

Similarity parameter for the cone boundary layer 

Index of mesh points in radial direction 

Jacobian of the transformation and/ or maximum 
number of mesh points in the radial direction 

Reference length 

Mass of a thermodynamic and/or momentum zone 

Mach number 

Radial momentum 

Axial momentum 

Distance normal to a surface 

Pressure 

Artificial viscosity function 
Radial coordinate 
Reynolds number 
Distance along a surface 

Coordinate velocity in the radial direction 
Safety factor for At c 

Safety factor for At d 
Temperature. 

Time 



At Times tep 

At c Sound-speed controlled timestep 

At^ Viscous-diffusion-controlled timestep 

U e Velocity component at outer edge of boundary layer 

U Radial velocity component 

U • Freestream velocity 

CO 

U Local particle speed (U = a. U 2 +W s ) 

W Axial component of velocity 

x Cartesian coordinate of abscissa in cross-sectional 

plane of calculation 

Cartesian coordinate of ordinate in cross-sectional 
plane of calculation 

Axial coordinate 

Freestream angle of attack with respect to the body 
axis 

Stream function 
Viscosity 

Kinematic viscosity (v=— ) 
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Dens ity 

Generalized radial coordinate 
Shear stress 

Generalized time coordinate 

Normal stress • component in the radial direction 

Normal stress component in the azimuthal direction 
* 

Normal stress component in the axial direction 
Velocity difference W-U 

CO 

Freestream condition 

Property -at outer edge of boundary layer 
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( ) n_1 

( > n - % 
( )* 


Reference condition for viscosity-temperature law 
Property associated with flow on a flat plate 
Property evaluated at t n 

II™ 1 

Property evaluated at t 

n-J- 

Property evaluated at t 2 
Non-dimensional property 


iv 



SECTION 1 


INTRODUCTION 

The purpose of this research effort is to develop a 
numerical method for computing steady, three-dimensional, 
viscous, compressible flow fields about aerodynamic bodies. 

The primary emphasis of this work is for supersonic axi- 
symmetric flow; however, the method should extend to three 
spatial dimensions and to the subsonic and transonic - flight 
regimes as well. - Axis ymme trie supersonic flow was -siudied^l' 
to facilitate development of the numerical method. The 
advantages of starting with axisymmetric flow are discussed 
later in this section. 

In 1969, under NASA AMES sponsorship, a research effort 

was initiated to apply the AFTON time- dependent two-dimensional 

computer code'*' to solve steady problems in three spatial dimen- 

2 3 

sions. The idea was based on the Equivalence Principle, ’ 
which states that for slender bodies at hypersonic speeds, the 
three-dimensional steady equations of motion for inviscid, 
flow reduce identically to unsteady equations in two spatial 
dimensions . 

The Equivalence Principle was extended in an ad hoc 
manner to a viscous flow through a model which permits viscous 
cross-flow, together with inviscid axial flow. Figure 1 
shows an ogive -cylinder body at angle of attack a with respect 
to the freestream flow direction. Leeward vortices are indi- 
cated in the figure. The axial coordinate z was made 
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( 1 ) 


proportional to a time- like variable t aecording to^ tbe 
relation 

z = U cos a t 

CO 

The two-dimensional Navier-Stokes equations were then solved in 
cross-sectional planes normal to the system's axis. These cross- ' 
sectional planes were moved at freestream speed., U^, from the 
leading edge to the trailing base of the body. The cross- . 
sectional flow at each time, t, corresponded to steady flow at 
the z-coordinate given by Equation (1) . It is noted that the 
.leeward vortices of Figure 1 have axes almost parallel to the 
z-axis; therefore, the cross-sectional planes contain the 
vortices and the vortex calculation should be under the scope 
of the method. 

Although all axial effects were neglected, this numerical 
procedure did calculate leeward vortices and produced other 
flow field results which were generally in accord with experi- 
ment «»' ■ -Numerical calculations were made and compared with 
experimental results for an ogive- cylinder and an airplane 
fuselage. Flow conditions were M^ = 1.98, = 4. 68x10 5 /ft 

and a - 10° for -the ogive -cylinder and M^ = 2.5, R^ -9.1x10 s /ft 
and a = 15° for the" fuselage configuration. Numerical- 
experimental comparisons were made of static pressure distribu- . 
tions on the body surface and contour plots of pitot pressure, 
total pressure, Mach number, and flow inclination at various 
cross-sectional planes normal to the system's axis. Maximum 
deviations of - numerical surface pressure's from corresponding 
experimental values were no more than 6 percent of free- 
stream dynamic pressure for both problems. Calculated pitot 
pressure contours were generally in accord with measured 
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contours on the windward side of the body. However, some 
differences were noted. Boundary layer separation and body 
vortex positions differed from experimental locations on the 
ogive -cylinder , and the shock induced by the fuselage canopy 
was predicted at a slightly different location. These differ- 
ences were considered attributable to neglect of axial viscous 
effects, exclusion of turbulence phenomena, and the approxi- 
mations introduced by the Equivalence Principle in describing 
the inviscid axial flow. Further details are presented in 
Reference 4. 

In this research effort, axial effects are incorporated 
into the AFTON numerical procedure and a numerical solution 
to the full Navier-Stokes equations is generated by iteration. 
The steady-axisymmetric equations of motion are solved using 
a radial one- dimensional time -dependent computer code. Compu- 
tational economy gained by one -dimensional calculation enabled 
the required experimentation with the finite difference tech^- 
nique to take place; hence, development of the complicated 
iteration method was made practical. The axisymmetric itera- 
tion method is applied to calculate the axisymmetric flow 
field about a cone-cylinder-flare body immersed in a supersonic 

airs tr earn. ' 

» ■> 

A discussion of the numerical method is presented in 
Section 2, results of the cone-cylinder-flare calculations are 
presented in Section 3, and the conclusions reached In this 
research effort are presented in Section 4. 
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SECTION 2 


THE NUMERICAL METHOD 

This section presents a general description of the 
iterative technique for solving the axisymmetric flow equations, 
formulation of the integral equations of motion for axisymmetric 
flow, and formulation of the finite difference approximations 
to them. 

2.1 - General Description of the Iteration Technique for 

Axisymmetric Flow 

Rather than apply the- Equivalence Principle in the ad hoc 
manner discussed in Section 1.0, it was decided to solve the 
full .Navier- Stokes equations for axisymmetric flow by iteration. 
The calculational procedure would be similar to that described 
in Section 1.0, except that the axial convective and axial 
viscous terms dropped in the Equivalence Principle formulation 
will be evaluated from numerical data of a previous iterate. 
Consider the axially symmetric body shown in Figure 2. The 
steady-state axisymmetric flow field will be made analogous 
to time- dependent one- dimensional flow in a plane containing 
an expanding cylinder. The cross-sectional plane of unsteady 
analogy moves downstream with free-stream velocity and the 
outline of the moving boundary is given by the trace of the 
original shape in the cross-sectional plane; the moving boundary 
has the velocity indicated in Figure 2. At station -z (see 
Figure 2) the steady- state flow field in this plane is analo- 
gous to the flow field about the expanding cylinder (of 
radius a) at the time, t = z/U^. Due to axial symmetry the 
time-dependent flow in each cross-sectional plane is radial. 

To make the analogy valid for a viscous fluid, no-slip flow 
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is imposed at the expanding cross-section. Far from the body 
free-stream conditions are imposed (i.e., E = E^, p = p , 

U- = 0 , W = V ) . 

CO 

The procedure for obtaining the first viscous iterate 

takes place in two steps. First, the inviscid solution or 

zeroth iterate is obtained. The flow field for the zeroth 

iterate is developed by dropping the viscous terms and axial 

convective terms from the equations of motion and applying 

2 3 

the hypersonic Equivalence Principle 3 directly. The first 
viscous iteration is obtained from the zeroth iterate by 
marching in z from the leading edge to the trailing base of 
the body and integrating the continuity, radial momentum, 
axial momentum and internal energy relations in a very special 
way. Axial convective and axial viscous terms are treated as 
source terms in the integration process and evaluated from 
numerical data from the zeroth iterate; linear interpolation 
in the inviscid field is used to evaluate the axial convective 
and axial viscous terms . All other terms are evaluated 
directly. In this way the mass, radial momentum, axial momentum, 
and internal energy are calculated throughout the axisymmetric 
flow field. After marching from the leading edge to the trail- 
ing base of the body, the first viscous iterate is obtained. 

The second viscous iterate is obtained by going back to 
the leading edge and recalculating the flow field by marching 
to the base of the body. Calculation of the second iterate 
requires evaluation of the axial convective and axial viscous 
terms from numerical data for the first iterate. The procedure 
is then repeated to determine successive iterations until 
convergence occurs. 

The iteration technique is further explained through the 
equations of motion presented in the next section. 
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2.2 .Formulation of the Finite Difference Equations 

for Axisymmetric Flow 

In this section finite difference equations are formulated 
for mass conservation, momentum conservation, and an internal 
energy theorem*. For illustrative purposes the finite differ- 
ence equation for mass conservation will he derived. Let us 
consider a fixed axisymmetric body immersed in an airstream 
with freestream speed U parallel to the axis of the body; the 
axial direction is denoted by z and the radial direction by r. 
The conservation of mass for time- independent motion in the 
Eulerian coordinates r, z is as follows: 


a(pW) _ 1 _ S(oUr) = n 

Sz r Sr 


( 2 ) 


where U is the local radial velocity, W is the' local axial 
velocity, and p is the local material density. We now trans- 
form Equation (2) to the unsteady (r,t) plane using the 
relations 



_ 1 _ _ 1 _ 
U St’ 


W - W + U 

CO 


(3) 


The transformed continuity equation is as follows: 


d(pr) S(pUr) = 

St 3r 


r at (pMr) 


(4) 


The left-hand side of Equation (4) represents the contin- 
uity equation for time- dependent one-dimensional radial flow. 

The right-hand side of Equation (4) represents a source term 
which accounts for the variation of the axial velocity W. This 
term must be considered known in the iteration process and will 
be evaluated from the flow field of the previous iterate in each 
successive iteration. 

The internal energy is defined as the difference between 
total and kinetic energies . 
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Equation (4) is formulated in the Eulerian coordinates 
(r,t); however., in the planar calculational process the trace of 
the body surface in the plane must distort with time. Hence, 
we are really interested in the continuity equation in a general- 
ized coordinate system |,t, where 

ts = T (5) 

r - f(?,r) (6) 

3f 

and f (§,o) — 5, while fg. - * ~ 1* 

S=0 T-0 


We now differentiate (5) with respect to t and r., respec- 
tively, with the results 


1 

0 



(7) 


Performing the same differentiation of Equation (6) yields 

f 

5g _ T 

St " £ e 



(B) 

(9) 


Consider the function 


SG 

Sr 

SG 

at 


G(|,t) whose derivatives are 


r 11 4- r — 
G S Sr + G t Sr 


Gj f e 
l § 


r . r 

G ? st + G t at 


-G_ -7T- + G 
? T 


( 10 ) 

( 11 ) 
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where use of Equations (7) - (9) produced Equations (10) and (11). 
Based on Equations (10) and (11) , the continuity Equation (4) 
transforms to the generalized coordinates (?,t) as follows: 

_ a(or) g f l_(pr )_ + a (p , Ur) . 

5? I dr ■ 35 

= " u“ {' If (pwr) s + (pwr) ^ (12) 

03 I ' 

where the grid velocity S is defined as 



Now the partial derivative — can be expressed as 

f 

as = _a_ /af \ = jfs 
a? ax \b%) h 


(13) 


(14) 


Based on Equation (14), the continuity equation becomes 


_a_ 

a? 


(orf |) 



pr(U-S) 


T if (» wrS ) 

CO 


jl_ _a_ 

U ar 



Now integrating Equation (15) over 


(15) 


_a_ 

a-r 



prf ^ d? + 


pr(U-S) : 


1 

U 



'J pwrf ? d| 


(16) 
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where Liebnitz’s rule is used to exchange differentiation and 
integration. The final step is to transform the ? integrals of 
Equation (16) to r integrals. We know that 


d| = Jdr 


(17) 


where 


M M 

3r Bt _ _1_ 

Bt Bt 
3r Bt 


Hence, the final integral continuity equation in generalized 
coordinates becomes 


B 

Bt 



-r 

-3 


i‘.f 

-2 


prdr 4- p(U-S)r 


+ p r(U m -W)S 

05 



pr (U -W)dr 

03 


ri 


(18) 


where Equation (3) is used to substitute for w in Equation (18) . • 

The corresponding integral equations for radial and axial 
momentum, and specific internal energy are presented in Appendix A. 
It is noted that the heat conduction term is left out of the 
specific internal energy equation of Appendix A. The Stokesian 
stress -rate-of-strain relationship, the equation of state for a ■ 
perfect gas, and a power law viscosity-temperature relationship* 
complete the system of equations. 


We now derive the finite difference analog to continuity 
Equation (18') . Consider a finite difference mesh along the 


*The viscosity- temperature relationship is presented in Section 3. 
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r axis with the grid positions r . , J.- The finite dif- 

ference mesh is illustrated in Figure 3. As seen from Figure 3 
the region of calculation is divided into two sets of zones, 
one for the calculation of thermodynamic variables such as den- 
sity, stress, internal energy, and the other for the calculation 
of momentum and velocity. Thermodynamic zones are indicated by 
solid lines and momentum zones by dashed lines. Thermodynamic 
properties are assumed homogeneous in solid zones and the momen- 
tum and velocity are assumed homogeneous in dashed zones. The 
AFTON 1 numerical procedure is employed to formulate the finite 

5 

difference equations for this continuum structure . Let the 

coordinate velocities S. and particle velocities (U. , WC) be 
' J J 

defined at the grid positions r^ and let the density Pj^ he 

defined at the zone centers r . tl =%(r.+r . ) . Based on the 

continuum structure of Figure 3, the finite difference approxi- 
mation to- Equation (18) for thermodynamic zone (r., r. ..) is' 

• 1 ' J J +-*- 


n 

m. 


n-1 
- m. 


.1 +% j +% A A 


-At 


= M + S 


(19) 


n-1 . 


where nu., .is the mass of zone (r., r. at time t n , nu,-. 

d+% J .1+1 J+i 

the zone mass at time ti ^ , At the times tep and 


is 


'A 

M = - 


n-1 

’j+l 




/ n-% c n- % \ n _ 

{ j+l ’ j+l “ P 


n-1 

j 


;U n-% _ s n-%\n-%; 
3 J / J 1 


_JL_ 

U 


n.-l 

v . t1 

J+l 


U - w n +\ s n + 

“ J+l/' J+l J+l 


p "' 1 |U_ -v i 


^2 1 1 


r . 
J 


s; 

J 


■(-20)- 
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+ 


Ml1-1 

- p j^ 


1 \/ ; r 1 ?,^+r^ / n-1 n-l> 

J - \ J ih+i- r i 


+ U -W' 

“ J+l 


n-1 -n-1^ 


n-1 n-l\ ] ) 

: j+r r j +%/J | 


( 21 ) 


rv-Tl r^n — "I <~n- 1 

The axial particle velocities W., W. , W. „ and densities 

1 J J+l 

~n <~n— 1 

Pj+i,’ p j+3- corres P on d to data from the previous iterate. The 

transport density ^ is defined as p 1 ? J} - = p^,i if the parameter 

J +1 J +i 1 +5 




-,n 


-+ 


n-1 


n-1 


U 


. S. ..f is positive and p..., = p., 3 if the parameter 

1+1 J+l 1+1 1+2 

n-1 


n ^ - S^.f is negative. The transport density p 


j+l j+l 


s 

n- 


1S 


j+l 

i. 


defined as p n ^ = p 1 ? , * if the parameter S 1 ? (U -W 1 . 1 ,?) is -positive 
s^ +1 *i+% r 1+1 v ~ 1+1' * 

and Pg 1 “ Pj + ^ if the parameter S^ + ^(U m -W^~^) is negative. 

j -j-p 

Considering the density to be homogeneous in -zone (r ^ , r_._^) 
yields 


n _n //n . n s , n n, 

p. pl “ 2m. ,, / (r. ,,+r.) (r.,--r.) 
l+% l+% 1+1 T 1+1 1 


( 22 ) 
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2.3 


Total Energy Conservation and Stability of the 
Finite Difference Equations 

The finite difference equations presented in Section 2.2 

and Appendix B are basically AFTON one -dimensional cylindrical 
5 

equations with known source terms. As in the case of all 
AFTON finite difference relations, the equations of mass con- 
servation, momentum conservation, the specific internal 
energy relation, and total energy conservation are redundant, 
i.e., any one of them can be deduced from the other three. 

More generally, any of the finite' difference equations for 
mass conservation, total energy conservation, and the inter- 
nal energy relation can be deduced from the other two taken 
together with the equation for momentum conservation. The only 
difference between the finite difference equations of Section 
2 . 2 and Appendix B and the usual AFTON equations is that the 
source terms must be carried along in deriving the total 
energy equation from the continuity equation, axial momen- 
tum equation, radial momentum equation, and the internal 
energy relation. 

Stability of the finite difference equations "of Section 
2.. 2 and Appendix B are' governed by the usual sound speed and 
diffusion' criteria; the source terms are known functions and 
therefore do not affect stability. The timestep for the cal- 
culation is governed by the smaller of the timesteps based 
on a sound speed criterion and a diffusion criterion. 

Ar . 

At = ■ _ m — ■■■ (23) 

C (U+c)S c 
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min 


2 v S , 
d 


( 24 ) 


where A is the sound- speed controlled timestep, At^ the 

viscous diffusion controlled timestep, Ar . the minimum radial 

— m:Ln 2 2 
increment, U the local particle speed (U = U 4W ) , c the 

local sound speed, v the kinematic viscosity, and S and S, 

c d 

are the safety factors for sound speed controlled and diffu- 
sion controlled timesteps, respectively. The parameters 
and were set equal to two. 
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SECTION 3 


RESULTS OF CALCULATIONS FOR A 
CONE-CYLINDER-FLARE BODY 

In this section a general description of the problem is 
given, the inviscid solution to the problem is discussed, vis- 
cous solutions to the problem are presented, and the effects 
of separation on the iteration procedure are investigated. 

1 {. 

3.1 General Description of Problem 

After a literature survey and consultation with the 
Project Officer, supersonic flow over a cone-cylinder-flare 
body was selected as the problem for numerical evaluation. 

A sketch of the cone-cylinder-flare geometry is presented in 
Figure 4. The cylinder diameter is D = 1.25 inches and the 
cone and flare have 10 degree half-angles. The freestream 

Mach number is 4.54 and the freestream Reynolds number (based 

on cylinder diameter) is R = 3.62x10^. This case was experi- 

7 50 

mentally studied by Kuehn at AMES Research Center. The 
quantities visually observed during the experiment were: 

(1) position of the separation point, and (2) the flow-field 
configuration, i.e., cone shock, flare shock wave-boundary 
layer interaction, and separated -.region at flare. The static 
pressure distribution was recorded along the cylinder -flare 
surface. These data will be compared with the numerical 
results . 

The reasons behind the selection of the aforementioned 
case are outlined as follows: 

(A) The equivalence principle is readily applicable to the 
slender sharp-cone configuration. 
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(B 1 ) Of the series of cone-cylinder-flare flows studied by 
Ruehn, the case selected had the highest Mach number, 
and therefore provides the simplest test of the applica- 
bility of the iteration technique. 

(C) The relatively short cylinder results in correspondingly 
few grid points in the streamwise direction, and there- 
fore minimizes computer time requirements . 

(D) The Reynolds number is high enough so that separation 
indeed occurs . 

(E) The Reynolds number is low enough so that transition 

probably takes place at some distance from the model. 

At a similar Mach number, transition is observed at the 

5 

flare at a Reynolds number of 5.35x10 , while there is 

, ' " “ ~ 5 * 

no sign of transition .at ! a Reynolds number of 4.97x10 , > 

Hence the entire boundary bayer for the selected case is 
laminar . 

(F) The flare angle of 10° is large enough to give an observ- 
able separated zone. 

(G) The boundary- layer thickness at the location where separa- 
tion is about to occur is approximately 0.042"', or 1/30 

of the cylinder diameter. Hence, a reasonably detailed 
picture of the flow field can be expected from the com-» 
putation . 

(H) At the station referenced in item (G) , the static pres- 
sure is about 90 percent of the free-stream pressure; 
there is no region in the flow field where the density 
is substantially less than the ambient density, hence, 
it will not be necessary to employ unusually high mesh 
point densities for the calculation. 

The medium considered was air with a po ly tropic Texpo- 
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nent Y of. 1.4 and with the following viscosity- temperature 

law , 

.98 


p, = p, 


./ E 


r \ E ' 


(25) 


5 2 8 

where u = 3.33x10 dyne- sec/cm , and E = 3.405x10 dyne-cm/ 
jt r 

gram* 


3.2- Mesh Boundary -Conditions-^- and Initial -Conditions Used 


The finite difference mesh for the cone-cylinder-flare 
problem was composed of 50 points and automatically distorted 
as the cylindrical trace of the body distorted in the cross- 
sectional plane. The -finite difference mesh in the vicinity 
of the body at axial station z = 8.29 era. is shown in Figure 
5. -The wedge- like radial region of calculation has an up- 
stream boundary at the radius r = 1.46 cm. and has its later- 
al boundary at-r = 7 cm. The zones are closely spaced in 
the vicinity of the body and coarsely spaced near the lateral 
boundary. About six zones exist in the boundary layer near 
the cone-cylinder junction; this is further discussed in 
Section 3.4.2. 


Boundary conditions along the radial wedge- like region 
of calculation were as follows: 

(1) No-slip flow on the upstream boundary. 

(2/).. i Freestream conditions on the lateral boundary, i.e., 
d — o ,E=E,U-0, W=U. 

CO 3 CO 3 9 CO 

(3) Axial symmetry on the side-walls of the wedge. 

The calculation did not commence from the leading edge 

of the cone, but started at a distance of z = .567 cm. down- 
* o 
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stream of the leading edge. At the cross-sectional plane 
corresponding to this downstream distance, impulsive initial 
conditions were imposed, i.e., freestream flow about the cone 
surface (p = p^ 3 E = E m , U = 0, W = W ) and no-slip flow at 
the surface (U = 0,.W = 0). 

3.3 Zeroth Iterate: The. Inviscid Flow Field 

The zeroth iterate, in an iterative computational pro- 
cedure to solve for the viscous, compressible flow about a 
cone-cylinder-flare configuration, is described in this sec- 
tion. The zeroth iterate consisted of the steady, inviscid 

flow field about this configuration subject to the assumptions 

2 3 

of the Hayes Equivalence Principle. J To solve this inviscid 

cone-cylinder-flare problem a hydrodynamic computer code, 

written in a generalized coordinate system, was developed 

for transient one-dimensional cylindrical flow, and the 

transient solution was converted to the steady inviscid 

flow field about the cone -cylinder- flare configuration via 

2 3 

the Hayes Equivalence Principle. 5 This computer code is 
designated "AFTON 1". The AFTON 1 computer code was devel- 
oped by programming the finite difference equations of 
Reference 5 to solve for flow fields on a finite difference 
mesh whose velocities can be prescribed in an arbitrary man- 
ner. The finite difference equations for the Lagrangian 
coordinate system were first programmed and debugged. While 
running the cone-cylinder-flare problem, it was found that 
the Lagrangian formulation produced zones so thin that the 
calculational timestep was reduced to a level which increased 
computer time requirements above the practical range. Hence, 
the generalized coordinate system was introduced and zones 
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were prescribed to be thicker as the problem progressed. 

Based on the AFTOH 1 computer code and the finite difference 
mesh of 50 points, the machine time required to solve the 
inviscid cone-cylinder-flare problem was 5 minutes on a Univac 
1108 computer. To attain this solution, whose results are 
discussed in the next two paragraphs, 2500 cycles (i.e., time- 
steps) of calculation were required. 

Results of the numerical calculations are compared to 

7 ' ' * ■ " 

the experimental data of Kuehn and with the known inviscid 
solution for the cone portion of the body. The numerical 
cone shock angle and cone surface pressure ratio (P/P ) were 
16.2° and 2.065, respectively, while the corresponding theo- 
retical inviscid values are 16.5° and 2.08. Further details 
of the comparison of the numerical results with experiment 
are -presented below. 

The numerical results for the inviscid cone-cylinder- 
flare body are compared to the experimental data of Kuehn in 
Figure 6. The body configuration is indicated in the figure 
(cross-hatched line) as well as the numerical cone shock locus 
flare shock locus, and surface pressure distribution. The 
experimental surface pressure distribution, cone shock locus, 
separation shock locus, and flare (or reattachment) shock 
locus are also indicated. It is seen from Fdgure 6 that the 
numerical shock angle is slightly smaller than the measured 
shock angle; this is in accord with physical reality because 
the boundary layer displacement effect is not present in the 
numerical inviscid calculations. The experimental pressure 
distribution is close to the numerical one upstream of the 
flare-corner interaction region and approaches the numerical 
pressures downstream of this interaction region; this result 
is physically reasonable since viscous effects are only im- 
portant in the neighborhood of the flare-corner. Finally, 
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the experimental flare (or reattachment) shock locus nearly 
coincides with the corresponding numerical shock locus . This 
result is at first somewhat disturbing, since the flow up- 
stream of the flare-corner has passed through the separation 
shock, and thus has a lower incident Mach number than the 
corresponding inviscid flow. Hence, the experimental flare 
shock angle should be greater than the inviscid angle. How- 
ever, due to separation the flow immediately downstream of the 
separation shock is flowing at a nonzero angle. Therefore the 
flare angle relative to this stream of fluid is less than the 
geometrical flare angle of 10° and this will produce a flare 
shock angle smaller than the inviscid one. The two counter- 
balancing factors can result in an experimental flare shock 
angle almost equal to the inviscid prediction. 

3.4 First, Second and Third Iterates: The Viscous Flow 

Field 

The finite difference equations for steady viscous, 
compressible, axisymmetric flow, which were formulated in 
Section 2.2, were incorporated into the AFTON 1 computer code. 
The revised AFTON 1 computer code was then capable of solving 
the viscous cone-cylinder- flare problem of Section 3. The 
solution was obtained through an iteration process based on 
the Equivalence Principle solution as the zeroth iterate. 

The revised AFTON 1 code was then checked out by performing 
hand calculations with the finite difference equations of 
.motion and comparing them with computer output for the cone- 
cylinder- flare problem. In addition to the hand calculations, 
the conservation properties of the AFTON 1 finite difference 
equations were checked. After seven cycles, i.e., timesteps, 
of computation, a zone-by-zone summation of the total energy 
(internal plus kinetic) of material in the region of 
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calculation was identical (within seven places) to the total 
energy that results from the conservation of total energy 
theorem for the entire region of calculation. 

Based on the AFTON 1 computer code and the finite dif- 
ference mesh of 50 points, the time required to solve for 
the first viscous iterate to the cone-cylinder-flare problem 
was 40 minutes on the Univac 1108 computer. To obtain this 
solution, whose results are discussed in the -remainder of 
this section, approximately 7000 cycles (i.e., timesteps) 
were required. This represents a considerable increase in 
computational time with respect to the zeroth iterate, which 
required 5 minutes and only 2500 cycles. However, a time- 
dependent solution to the unsteady Navier- Stokes equations 
for axisymmetric flow would require about 12 hours on a 
Univac 1108 computer, so that the inherent economy of the 
procedure is clearly indicated through the first two iterates. 

3.4.1 Flow Field and Pressure Distribution 

A velocity vector plot of the cone-cylinder-flare flow 
•field after three viscous iterations is shown in Figure 7. 

Each column of vectors represents the -solution in one cross- 
sectional plane. The vectors in each column emanate from 
the mesh points along the radial line in that cross-sectional 
planed, their magnitudes are proportional to the local particle 
speeds. Curve 1 denotes the calculated cone shock locus and 
Curve 2 denotes the calculated flare shock locus. Shock trans- 
itions were calculated with an artifical viscosity, i.e., 

Q- function, proportional to the square of the velocity diver- 
gence. However, the Q- function was limited only to the 
neighborhood of the shock transitions. A thin boundary layer 
is seen as a dark "line near the body surface; the calculated 
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boundary layer thickness was approximately 2% of the cylinder 
diameter near the cone-cylinder junction. 

The numerical results of the zeroth and first iterates 
are compared to the experimental data of Kuehn in Figure 8. 

The body configuration is indicated in the figure as well as 
numerical cone .shock loci, flare shock loci, and surface pres- 
sure distributions for both the zeroth and first iterates. 

The experimental surface pressure distribution, cone shock 
locus, separation shock locus, and flare (or reattachment) 
shock locus are also indicated. The experimental cone shock 
angle is 17.6°; by comparison the first iterate cone shock 
angle is 17.0°, while the cone shock angle of the zeroth 
iterate is 16.5°. This represents a 3.4% error in cone shock 
angle for the first iterate and an 8.5% error for the second 
iterate. The flare shock angle has remained almost constant 
in going from the zeroth to the first iterate. 

The surface pressure distributions of the zeroth and 
first iterates are almost identical and clearly do not match 
the data in the region of separation at the cylinder- flare 
junction because the boundary layer computed in the first 
iterate did not separate at the cylinder- flare junction. 
Subsequent iterations would be required to determine the sub- 
sonic recirculation 'region associated with -separation. 

3.4.2 The Cone- Boundary Layer for Infinite! Prarid.tl Number 

As was mentioned in Section 2.2, heat conduction was 
left out of the internal energy relation solved in this nu- 
merical investigation. Hence, the numerical solutions gen- 
erated correspond to infinite Prandtl number. Unfortunately, 
to the best of our knowledge, there are no available boundary 
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layer theory predictions to compare to for flow along a cone 

at infinite Prandtl number. However, the Chapman-Rubesin flat 

9 

plate boundary layer velocity profile, corrected for axial 

10 

symmetry by the Mangier transformation, is independent of 
the energy equation if plotted as a function of the similarity 
parameter 


f 


/up U Cs 
A e e e 


(26) 


where ij> is the stream function for axially symmetric flow, 

U i-s the velocity at the boundary layer edge, p is the 
density at the boundary layer edge., is the viscosity -at 
the boundary layer edge, C is the Chapman-Rubesin constant, 
and s is distance along the cone surface. The Chapman- 
Rubesin boundary layer velocity profile for flow along a : _~ 
cone is discussed further in Appendix C. 

The numerical results for the cone boundary layer 
velocity profile, corresponding to the first and second 
iterates are compared to the Chapman-Rubesin profile in 
Figure 9. The numerical data correspond to an axial station 
z of 8.73 cm. and a Chapman-Rubesin constant C of i.8875. 

The velocity profile for the first iterate has an inflection 
point and has a maximum deviation from the theoretical pro- 
file of no more than 17%. However, the boundary layer 
profile for the second iterate coincides with the Chapman- 
Rubesin 'prediction. Hence, after two viscous iterations 
the cone boundary layer has been calculated .correctly. 
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3.5 


Separation at the Cylinder-Flare Junction 

One of the principal objectives of this research effort 

is to investigate the feasibility of predicting flow separation 

7 

through the iteration procedure. Experimentally, it has been 
shown that the flow separates at the cylinder-flare junction 
and an axisymmetric vortex forms there. In this section the 
numerical results of five viscous iterations are examined in 
the neighborhood of the cylinder- flare junction to determine 
the extent to which the method can predict separation. 

Development of the flow field in the vicinity of the 
cylinder-flare junction is presented in the form of velocity 
vector plots in Figures 10, 11 and 12. Figure 10 presents 
the flow 'field after three v-iscous iterations, Figure 11 
after four viscous iterations, and finally Figure 12, after 
five viscous iterations. The + symbols in each figure repre- 
sent the body surface. It is seen from Figure 10 that after 
three iterations the flow is only slightly perturbed at the 
cylinder-flare junction, and definitely is not separated. 

Figure 11, after four viscous iterations, shows the - flow to 
be near separating but not quite separated. Finally, Figure 
12, after five viscous iterations, shows that separation has 
been initiated. 

Based on the above results, it seems clear that the 
calculation of flow separation through successive iterations 
is possible but may not be -practical. It is anticipated that 
approximately fifty ■ iterations may be required to completely 
define the vortex at the cylinder- flare junction. 

The cone-cylinder-flare problem could not be continued 
past the fifth viscous iteration. During computation of the 
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sixth viscous iterate, large oscillations were recorded along 
the cone surface, which finally 'produced a negative mass-- 
i.e., a "numerical singular point." It is believed that dis- 
cretization error, defined as the difference between the 
exact value of a property of the continuum and the value de- 
termined from a discrete approximation to the continuum, is 
the '.primary cause of termination of the iteration procedure. 
Past experience has shown that discretization errors produce 
small spatial and timewise oscillations in the numerical flow 
field. Flow field oscillations in the cone-cylinder-flare 
problem 'built up fluctuations in ’the source terms of the equa 
tions of motion in each successive iteration. The fluctuat- 
ing source terms eventually caused the failure on the cone 
surface during the sixth viscous iteration. It is believed 
that the iteration procedure could be continued if the dis- 
cretization oscillations were -removed from the flow field. 
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SECTION 4 


CONCLUSIONS 

Based on the results of Section 3, two principal conclu- 
sions emerge regarding the iteration technique , -for solving the 
full Navier- Stokes equations. 

(1) The iteration technique is capable of correctly pre- 
dicting details of the flow field which have little 
or no upstream influence, such as a boundary layer 
or a shock wave. 

(2) Although theoretically possible, it is impractical to 
calculate separations which lie outside the plane of 
calculation: through the iteration procedure. The 
cylinder- flare junction vortex has its axis parallel 
to the cross-sectional planes of calculation; hence, 
it lies outside these planes. 

The above conclusions indicate that the application of 
one-dimensional time-dependent methods to solve for steady 
plane two-dimensional and/or axially symmetric flow fields 
are of limited value. These one- dimensional methods will be 
capable of predicting only those flow phenomena with little 
or no upstream influence, i.e., a boundary layer or a shock 
wave. Vortices, in plane two-dimensional and axisymmetric 
flows will always lie outside the planes of calculation; 
therefore, they are beyond the scope of one-dimensional 
methods . 

Conclusions (1) and (2) also clearly indicate that the 
application of two-dimensional time- dependent methods to solve 
for steady flow fields in three spatial dimensions is of great 
value. In a three dimensional flow problem, one can always 
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select a two-dimensional plane of calculation which will be 
nearly normal to the. axes of vortices in the flow field. 

There are two excellent examples of three-dimensional flows 
which are definitely under the scope of this iteration tech- 
nique. The first is a body at angle- of- attack and the second 
is a three-dimensional wing. 

The body at angle- of-attack problem is a simple one 

for the iteration technique. Experimentally determined axial 

1 JL . o 

velocity profiles on the lee side of a 10 half-angle cone, 
at 12° angle of attack, freestream Mach number 7.87, and Rey- 
nolds number 8.51x10^ (referred to the cone base diameter), 
and boundary- layer -like. The recirculation region occurs in 
cross-sectional planes normal to the body's axis. Applica- 
tion of the iteration technique in cross-sectional planes 
will surely be successful. 

The three-dimensional sweptback wing problem is more 
difficult for the iteration technique, but definitely feasible. 
This problem generally has two types of vortices. In a plane 
normal to the wing span, i.e., a cross-sectional plane, the 
principal vortex will exist. Furthermore, due to the presence 
of the wing tip, a secondary vortex could exist in a plane nor- 
mal to the free-stream flow direction. To solve the wing prob- 
lem the integration must be done in two different directions 
in the iteration procedure. First, two iterations are required 
in cross-sectional planes of the wing. These iterations will 
define the principal vortex. Secondly, a third iteration is 
required in planes normal to the free-stream flow direction 
to determine the secondary (tip) vortex and other wing tip 
effects. It is believed that in three iterations, as described 
above, the flow field about a three-dimensional swept-back 
wing can be calculated. 


26 



APPENDIX A 


INTEGRAL EQUATIONS FOR AXISYMMETRIC FLOW 


In this appendix the integral equations for momentum 
conservation and the internal energy theorem are presented 
for axisymmetric flow. These relations are first presented 
in differential form for an Eulerian coordinate system and 
then shown in their integral form for a generalized coordinate 
system. 

Consider a fixed axisymmetric body immersed in an air- 
stream with freestream speed parallel to the axis of the 
body; the axial direction is denoted by z and the radial 
direction by r. The radial momentum equation-, axial momentum 
equation, and internal energy relation in Eulerian cylindri- 
cal coordinates r and z are as follows : 
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where E is the specific internal energy, U is the local radial 

velocity, W is the local axial velocity, and p is, the local 

material density. The parameters o' , t , cr and o are 

rr zr 0 0 zz 

the stress components in the radial, tangential, azimuthal, 
and axial directions, respectively. Equations (Al) , (A2). and 
(A3) are then transformed to the unsteady (r,t) plane using 
the relations 
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The transformed equations are formulated in the 
Eulerian coordinates (r,t); however, in the planar calcula- 
tional process the trace of the body surface in the plane 
distorts with time. Hence, we are really • interested in the 
equations of motion in a generalized coordinate system 5, t, 
where 


t — t (A5) 

r = f-(?,T) (A6) 

and f(§,0) = |, while f ^ = ('-— ) t= q = 1 • Proceeding 

through the same derivation as outlined in Section 2.2 for 
the continuity equation yields the remaining integral rela- 
tions . 
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r 2 . 

where the bracket [ ] implies that the term contained is evalu- 
ated at the limits r^ and.^* Equations (A£)j- (A8i) , and (A9) 
in conjunction with Equation (18) represent the four conservation 
relations solved numerically for the dependent variables p, E, 

U and W. 
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APPENDIX B 


FINITE DIFFERENCE EQUATIONS FOR 
AXISYMMETRIC FLOW 


The finite difference analogs to integral momentum 
Equations (A5) and (A6) and the integral internal energy 
Equation (A7) are presented in this appendix. Consider the 
finite difference mesh along the r axis with the grid positions 
r j , j=l, J. The finite difference mesh is illustrated in 

Figure 3. Let the coordinate velocities and particle velocities 
(U^. ,Wj ) be defined at the grid positions r ^ and let the density 
and energy E^ be defined at the zone centers r j + p = 
^(rj+r. + ^. Based on the continuum structure of Figure 3, the 
finite difference approximations to the momentum equations 
(Equations (A7,) and (A8)) are written for the momentum zone 

(r. i, r j), while the finite difference approximation to the 
3 - i> 3+-^ y 

• internal energy relation (Equation A9,)) is written for the ther- 
modynamic zone (r_. , r_. + ^,). These relations are as follows: 
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where T n is the shear stress of zone r., r.,, evaluated at 
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time t.. 2 and t to the shear stress of zone (r.-r..,) at, 
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derivatives of z in the formulae for T n and T n ^ . The 

zr . , j zr 1 

n -l l+% J + ^n 

stress t 2 is taken as the arithmetic average of t 

zr, xl ° zr -.i 

n _l j+% 1+^ 

2r j+%' 


and 


Axial Momentum 


,,n , Ul” 1 

M -M 

Z A Z A A A A A 

— ® Mz + Sz + Fz + Gz 

At 


(B6) 


where M n is the axial momentum of monentum zone (r. ls r. j) 
Z J! l~i> J+'S 

J ^ /\ /\ A A 

at time t , and Mz, Sz, Fz, and Gz are defined as follows: 


A 

M = - 

- z 




1 W 

n-%,"i+l i v n-^, 

•W , ) ' “j-% { 


n-% / J J 1 > 


_1 

U 


, n-% , J . 

,1^.1 ( 2 > 

=0 (I J+% 


- M- 


n-% , J 
s . , '■ 




) (B7) 


A 


. i 

YT7- { (u -w") Sf? 

U . At ), “J1 


n 


n 


~n f 

P j+^ 


r j+% +r j^,n n 


\ f LI tlx 

)(r j+%' r j ) 


n , n 
r • +r • i 

. ~ n / J J~t \ / T . n „ r n \ 

+ W 2 )(r j j-%% 


- (u -w?' 1 ) w n_1 

05 1 J 


n-1 n-1 

!'~n-l r r j+% +r j w n-1 n-1, 

P j+% ( 2 >< r j+ % J 7 


n- 1 , n- 1 

- r • +r • v n t 
, ~n-l, J J — ”2 N / n-i jn-l x 

■ + W — —) ^ 


l 


W) 


34 



r r 


A 1 
Gz = 1 


U At 

OO ,1 


.n- 

sr 

Z’Z , 
•]!% 


n n n n 

r ■ + j,+r ^ r • v r . i 

f I 2 J % „fl-\ r ti' / j j — *2\ / x\ n N , 

( 2 )(r. +% -r j) + 0z . (— ^ )(r r r % )J 


.n-1 

a/ 

22 - , i 

J+% 


( r n_ ^+ r n_ ^W r n ''I_ r n '' 4 - -nE 1 "^ '(. 

C j+% +r j } (r j+% r j } + - ( 


n-1, n-1 

r j +r i-4 


J _; 2 n , n-1 n-1 


zz . 


-or - n-% r i. „ 1+1 V. n-%. n-% r J i-U j 
j j+% ( '2 ' } zz — ^ o ^ 


1 ) n-% 




l J-% 


2 

-lr _n-Ac 

2 H 


) ( J 


u is the axial normal stress of zone (r.,r.,,) at 
zz.,, v j ’ j-HL 


where a 

zz. ,1 
n ^ 

time t and <7 is the axial normal stress of zone (r ls r.,,) 

zz. ,i v .1 ’ j+1 

1 j 4% J j 

at time t . The stress a is evaluated from the sum- of 

zz. (1 

n J + "2 n _ j. 

a pressure term at t and a viscous term at t 2 ; velocity data 
from the previous iterate are used to determine derivatives of 


z in the a 


n formula. The stress o- n 2 


ZZ « . 1 

3+h 

the arithmetic average of a 


n 

ZZ * . 1 

J+% 


and ct 


zz. 

n-1 


is evaluated from 


+% 


zz. , 

J+% 


~ r , P. 

J 2 




Internal Energy Relation 


n n n-1 n-1 

i +% i +% i i +% A A A A ,A 

— i ~ 3 — J 2 — J 2 — J— 3 a = M„ + + W_ + G + 

At E E E Ez ER 


(BIO) 


where E 1 ? is the specific internal energy of thermodynamic 
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where the transport energies E. . and E are defined in a 
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of Section 2.2. 
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The specific internal energies - , E? \ correspond to data 

J +2 J+t> ’ 

from, the previous iterate. 
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where the particle velocities ^j+l’ ^ correspond 

to data from the previous iterate. 
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APPENDIX C 


CHAPMAN-RUBES IN BOUNDARY LAYER 
VELOCITY PROFILE ON A CONE 

The Chapman ^-Rubes in flat plate boundary layer velocity pro- 
file,^ corrected for axial symmetry by the Mangier transformation, 10 
is independent of the energy equation if plotted as a function of the 
similarity parameter 


(Cl) 
s 

where is the stream function for axially symmetric flow, U^ is the 
velocity at the boundary layer edge, is the density at the bound- 
ary layer edge, p, is the viscosity at the boundary layer edge, C is 
the Chapman -Rube sin constant and s is distance along the cone surface. 
In this appendix the Chapman-Rubes in boundary layer velocity profile 
for a cone is derived. This derivation starts with a flat plate geo- 
metry and then the results are transformed to a cone. 

The Chapman-Rubesin velocity profile for a flat plate is based 
on the principal assumption that the viscosity in the boundary layer 
is proportional to the temperature in the boundary layer. The follow- 
ing relation is assumed for the viscosity-temperature relationship: 


f = 



'p; CUs C 
e e 



T 

C ^T~ - 
e 


(G2) , 


where p, is the local viscosity, T the local, temperature and T^ the 
temperature at the boundary layer edge. 
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Chapman and Rubesin then introduce non-dimensional forms 
of velocity, temperature, stream function, density, and so on, 
as follows: 


U" 


_U 

Viu 5 

ee 





£ • < -*■ 


p 'U L C 
e e e 



(C3) 


where L is a characteristic distance, n is the distance normal 

i 

to the surface and the subscript "p" denotes a property asso- 
ciated with the flat-plate flow. The absence of a subscript 

i 

"p" denotes a property which is invariant to the geometry of 
the system, i.e., the same for a flat plate and a cone. In 
terms of the non-dimensional variables of Equations (C3), the 
boundary layer momentum equation becomes 






(C4) 


where the usual independent variables of non-dimensional 

distance along the plate s and non-dimensional .distance 

normal to the plate n have been replaced by the non-dimensional 

it P 

stream function i|r and non-dimensional distance along the plate 


Sp according to the von Mises transformation. 


12 


Equation (C4) 

does not have the temperature or density in it; hence, it is 
independent of the energy equation. 


Partial differential equation (C4) is then reduced to 
a total differential equation through the similarity parameter 
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(C5) 



and the assumption that 

U* = U*(f p ) (C6) 


Introduction of Equations (C5) and (C6) into Equation (C4) 
yields 


f p dU‘ d dU*x 

2 ir = dr (u dr } 

p p P 


(C7 ) 


Let us now assume that is a function of another dummy 
variable 11 and that the 'non-dimensional velocity profile is 
given by 


u 



(C8) 


Introduction of Equation (C8) into Equation (C7) yields the 
final differential equation 

I M IT 

2f +ff =0 (C9) 

p p p v 

Equation (C9) is identical to the Blasius flat plate equation 

13 

for incompressible flow. Therefore, according to Equation 

(C6) and (C8), the Blasius solution for the non-dimensional 
velocity profile on a flat plate represents the compressible 
profile if plotted as a function of the Blasius stream function 
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It is now necessary to transform the known flat plate 

velocity profile to suit flow along a cone. For this purpose 

the Mangier transformation"^ is used. The flat plate stream 

function f and distance' along the plate s transform to cor- 
P P 

responding cone properties according to the following relations: 


* 


s 





(CIO) 


where 9 c is the cone half-angle. The Blasius stream function 
f then transforms to the cone geometry as follows : 


f s sin 6 



(Cll) 


where f is defined, given by Equation (Cl). Hence, if the 
Blasius velocity profile is plotted as a function of f (given 
by Equation (Cll)), the boundary- layer velocity profile for 
compressible flow on a cone results. 
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Figure 1. Flow Field About an. Axisymmetric Body at Angle of Attack. 
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EQUIVALENT 

THREE DIMENSIONAL STEADY PROBLEM ' TWO DIMENSIONAL UNSTEADY PROBLEM 


Figpre 2. Application of the Iteration Technique to an Axisymraetric Body. 
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Figure 3. Schematic of Three Thermodynamic Zones and Two Momentum Zones for Radial One- 
Dimensional Time-Dependent Motion; solid lines indicate thermodynamic zones; 
dashed lines indicate momentum zones; a indicates cylinder radius. 
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Figure -4. 


Schematic of Cone-Cylinder-Flare Geometry Numerically 
Investigated. 
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CARTESIAN COORDINATE, y. cm 


Trace of cone in cross-sectional 
plane of calculation 


a - 1.46 cm 
A9 = 1° 


CARTESIAN COORDINATE, x, cm 


Finite Difference Mesh Along a Wedge-Like Radial Region of Calculation in the (x,y) Cross-section 
Plane 'at an Axial Station z of 8.29 cm; a denotes the radius of the trace of the body in the 
cross-sectional plane; A0 corresponds to the azimuthal angle of the wedge-like radial region; 
the flow is symmetrical with respect to 9 but varies in the radial direction. 
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V Experiment Cone Shock 

□ Experimental Separation Shock 

^ Experimental Flare (or reattachment) 
Shock 



o o o 


Axial Distance, z, cm 


Figure 6, 


Comparison of Inviscid Numerical Results with Experimental Surface Pressure Distributions 
and Shock Loci; M =4.42, R = 1.14xlCr per cm; solid lines indicate numerical results, 

Clj r Oft *■ * * 


cross-hatched line indicates body geometry, and symbols indicate data points. 
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Surface to Free Stream Pressure Ratio Radial Distance, r, cm 


0 Experimental Cone Shock 
.□ Experimental Separation Shock 
^ Experimental Flare (or reattachmeht) 



"• ”t- - r Axial Distance, z, cm ' ■ 

iFigure 8. Comparison of Numerical and Experimental Surface Pressure Distributions and Shock 
loci: M =4.42, R = 1. 14xlCr/cm; solid lines indicate numerical results for first 
' viscous iterate; dashed lines represent inviscid solution or zeroth iterate; cross 
hatched line indicates body geometry. 
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Figure 9. 


Comparison of Numerically Calculated Boundary Layer Velocity Profiles on the 
Cone Surface with Boundary Layer Theory. 
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Figure 10, Velocity Vector Plot of the Flow Field in the Neighborhood of the Cylinder- 
Flare Junction; vectors correspond to numerical data from the third iterate 
+ symbols indicate the cylinder-flare geometry. 
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